Discrete Nonlinear Schrodinger Equations with arbitrarily high order nonlinearities 
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A class of discrete nonlinear Schrodinger equations with arbitrarily high order nonlinearities is 
introduced. These equations are derived from the same Hamiltonian using different Poisson brack- 
ets and include as particular cases the saturable discrete nonlinear Schrodinger equation and the 
Ablowitz-Ladik equation. As a common property, these equations possess three kinds of exact ana- 
lytical stationary solutions for which the Peierls-Nabarro barrier is zero. Several properties of these 
solutions, including stability, discrete breathers and moving solutions, are investigated. 
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I. INTRODUCTION 

The discrete nonlinear Schrodinger (DNLS) equation 
appears ubiquitously 1] throughout modern science since 
it represents one of the simplest equations in which the 
combination of dispersive effects with a cubic nonlinear- 
ity leads to localized solutions of soliton type. Most no- 
table is the role it plays in understanding the propagation 
of electromagnetic waves in glass fibers and other optical 
waveguides 0. More recently the DNLS has been used 
as a tight bindingmodel for Bose-Einstein condensates 
in optical lattices 0,0. From a physical point of view, it 
is of interest to study the effects of including high order 
nonlinear terms (higher than cubic) in the equation on 
discrete solitons. These terms appear in different physical 
contexts such as Bose gases with hard core interactions 
in the Tonks- Girardeau regime Q and low-dimensional 
Bose-Einstein condensates in which quintic nonlinearities 
in the nonlinear Schrodinger (NLS) equation are used to 
model three-body interactions [||. A self- focusing cubic- 
quintic NLS equation is also used in nonlinear optics as 
a model for photonic crystals Q. 

For the continuous NLS equation with attractive in- 
teraction it is well known that the higher order nonlin- 
earities (higher than cubic) lead to the collapse in a fi- 
nite time (blow up) if the norm exceeds a critical value, 
even in the one-dimensional case. The interplay between 
dimensionality and the order of nonlinearity has indeed 
been used in the past as a way to investigate collapse in 
low dimensional nonlinear systems 8] . Although in a dis- 
crete NLS system true collapse cannot occur, due to the 
conservation of the norm, it may be possible that some 
of the features observed in the continuous NLS system 
about localized solutions may also exist at the discrete 
level. In particular, it is known that in the ID continuous 
NLS equation with high order nonlinearity (for example 
quintic) there exists only one localized solution for each 



value of the norm (critical norm), the so called Townes 
soliton Jgj, which separates collapsing and decaying so- 
lutions while being marginally stable against decay or 
collapse. In the presence of an external field, for exam- 
ple a periodic potential, it is possible to stabilize such 
solutions of continuous NLS with higher order nonlin- 
earities against decay, extending the existence range of 
localized solutions from a single value of the norm to a 
whole interval. Since the discrete NLS can be viewed as 
a tight binding model of the continuous NLS with a peri- 
odic potential, it is of interest to investigate the existence 
of discrete, stable localized solutions when higher order 
nonlinearity are introduced in the DNLS. 

The aim of this paper is twofold. Firstly, we introduce 
a general DNLS equation with arbitrary higher order 
nonlinearities which in the appropriate limits reduces to 
the integrable Ablowitz-Ladik (AL) equation Hi| and to 
the cubic DNLS with a saturable nonlinearity |l2| . These 
two equations have the remarkable property that they 
possess different analytical stationary solutions, both pe- 
riodic and localized (solitonic), which are exact to all 
orders of nonlinearity. These solutions exist for specific 
values of frequency and nonlinear parameters and have 
been shown to be stable under small perturbations. Sec- 
ondly, we investigate the effects of higher order nonlin- 
earity on the stability and mobility of localized solutions 
of a generalized DNLS. In particular, we compute the 
Peierls-Nabarro (PN) barrier and perform direct numeri- 
cal integration to show the existence of moving solutions. 

The plan of the paper is as follows. In Sec. II we first 
show that the AL model and the saturable DNLS 
model [12( can be obtained from the same Hamiltonian. 
In Sec. Ill we extend these ideas to Hamiltonians with 
arbitrary nonlinearity and obtain several higher order 
DNLS models. We then obtain a number of (time) sta- 
tionary, spatially periodic as well as localized solutions. 
In Sec. IV we study the various properties of these solu- 
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tions. In Sec. V we examine the question of the existence 
of moving solutions in DNLS models with higher order 
nonlinearities. Finally in Sec. VI we point out the possi- 
ble implications of our results and some open problems. 



II. THE MODEL 

In a classic paper, Ablowitz and Ladik showed that 
one of the variants of the DNLS equation given by 

ty„ + (1 + \^ n \ 2 )[ip n +i + tpn-i] ~ 2V> n = 0, (1) 

is integrable. In Ref. ^lj a model was proposed which in 
one limit goes over to the DNLS model while in the other 
limit it goes over to the integrable AL model. Recently, 
we were able to obtain exact periodic solutions of the 
DNLS equation with a saturable nonlinearity ^3] 

v\ib I 2 

tyn + (lpn+1 + i>n-l - 2V>n) + : , " = 0, (2) 

1 + H\ip n \ 

which is an established model for optical pulse propaga- 
tion in various doped fibers 0. In Eq. (0, t/>„ is a 
complex valued "wave function" at site n, while v and /i 
are real parameters. 

We point out that the two equations, i.e. Eqs. JU and 
J2J can both be derived from the same Hamiltonian Ti. 
given by: 



A' 



will assume fi ^ and perform the transformation 
■s/Jlipn —> ipn- This will replace \i by 1 and ^ by u, thus 
rendering the problem a one parameter problem. 

For the A = case the equation of motion, Eq. 
(0), can be written as 

*(1 + \*p n \ 2 )lp n + (1 + \lp n \ 2 )(i> n+ l + tp n -l ~ 2lp n ) 



+ I>|V»|V» = > (7) 

while for the A = 1(= fx) Eq. JHJ becomes: 

# n + (l+|V'n| 2 )(V'n+l+V'n-l-2?/'r l ) + ^l'/'™| 2 V'n = 0. (8) 

Note that v = 2 in Eq. (JSJ gives the AL equation [To| . 
Eq. (JTJ. 

In both cases a conserved power V can be written: 



N 1 

P = V^_ln(l+A|W 2 ), A->0 or A = /i = 1 . (9) 

n— 1 

The difference between the two cases is the presence of 
i|i/> n | 2 in the factor in front of the time derivative term 
in Eq. 0. However, for a stationary solution [i.e. only 
exp(— iujt) time dependence: iip n — unpn] Eq. (JJJ) can be 
written as 



|V»n - V'n+ll 2 \4>v 



ilp n + (1 + |^„| 2 )(^„+i + Ipr, 



2^n) 



and the equation of motion in both cases is 



(3) 



(4) 



The difference in the equations of motion comes from 
a different definition of the Poisson bracket (PB) and 
consequently a different definition of the time derivative. 
The Poisson bracket structure in both the cases can be 
compactly written as 



N 



dU dV dU dV 



' t dipn dip* dip n d^Pr 



)(l + A|V n | 



(5) 



On using Eqs. © and Eq. |gj for A = Nj^ and itp* 
are conjugate variables) it yields Eq. (J2J) [12] through 
Eq. J2J, while A = /j (ip n and are non-co njug ate 
variables) yields the equation introduced in Ref. 11] 

ilpn + (1 + fl\lp n \ 2 )(lp n+ i +1pn-l ~ tyn) + u\lp n \ 2 1p n = 0. 

(6) 

Notice that for fj, — > both Eq. @ and Eq. iJBJ reduce 
to the ordinary DNLS. Therefore, in the following we 



+ (^ + w)|V„| 2 V™ = 0. 
which is identical in form to Eq. (JSj. 



(10) 



The exact solutions to Eq. J7J given in [l^ all are 
stationary solutions "rotating" with the frequency 



uj = 2 — v . 



(11) 



Inserting this frequency into Eq. (|1U|) gives the AL equa- 
tion: 

ilpn + (1 + |Vn| 2 )0Wl + l/fn-1 ~ tyn) + ^nflpn = . 

From this is clear that exact stationary solutions of the 
saturable DNLS equation ^| are also stationary solu- 
tions of the AL equation. Analytical expressions for sta- 
tionary and moving solutions of the AL equation were 
given by Scharf and Bishop in Ref. (note that their u> 
differs from ours by a factor of two, see Eq. (fill) . We also 
remark that, because of the frequency relation lu + v = 2, 
Eq. © has stationary solutions of the AL only if v = 2, 
i.e. when it reduces to the AL equation. In the following 
we shall demonstrate the existence of analytical station- 
ary solutions of the AL type also for generalizations of 
Eqs. P)l. with arbitrarily high order nonlinearities. 
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A. High order nonlinearities 



III. EXACT ANALYTICAL SOLUTIONS 



We now generalize the nonlinear part of the Hamiltonian 
H in Eq. © by replacing ln(l + |-0n| 2 ) by: 



ln(l + /(|Vg 2 )), 



(13) 



where the function f(x), a polynomial of degree p+l, is 
given by: 



/(*) = £■ 



-3+1 



(14) 



3=0 



Here ao is always 1. Having already considered the case 
// — > 0, and having chosen fi = 1, we can next choose two 
different values of i/ in Eq. J2J. Therefore, instead of Eq. 
© , the Hamiltonian TL is now given by: 



N 



+ ^ln(l + /(|^| 2 )) 

and the generalized Poisson bracket by: 
[U,V] = 



(15) 



N 

E 



[1 + A/(|V„| 2 )], (16) 



with A = or 1. Note that the equation of motion is still 
given by Eq. Q. 

We notice for the A = case that a transformation 



-iAut 



i\) n would add Auj to the coefficient of ip n in 



the equation of motion and subtract the same from the 
coefficient of the \ip n \ 2 term in the Hamiltonian. There- 
fore, a v' different from v will only shift the frequency. 
We remark that some of the recently discussed models 
[TH fall in this category and are essentially equivalent 
to the saturated DNLS model [l^. Thus, both in this 
section and in Sec. IIIA the effect will only be to change 
the frequency u> by Auj. This will, however, not be the 
case for A = 1 as discussed here and in Sec. IIIB. 

We also note that with higher order nonlinearities too, 
besides the Hamiltonian Ti., the power 



N 



V = J2\^n\ 2 , A = 0. 



or 



N 



P = ^ln(l + /(|^ n | 2 )), A = l 



(17) 



(18) 



n=l 

is a conserved quantity. 



The main objective here is to find stationary solutions of 
the type obtained in Ref. [l2T | but with an additional am- 
plitude factor A in the equation of motion derived from 
the generalized Hamiltonian H given by Eq. I|15|) . I n 
particular, as in Ref. 0] we try to obtain three different 
types of solutions. The type I solution is given by 

, = A sn(M e - iM+S)dn(/3(n + c)) m) (19) 
cn(p, m) 

where the frequency u> is given by Eq. (flip while the two 
equations determining m and (3 are 



2K(m) dn(/3, m) 

P = — — — , v — 2- 



N„ 



cn 2 (/3, m) 



(20) 



Here N p denotes the number of sites in one period of the 
system. On the other hand, the type II solution is given 
by 

^ = AV^^f\e-^+ s )cn{l3[n + c),m), (21) 
dn(p, m) 

with the same frequency ui as given by Eq. i|ll|) and the 
two equations determining m and (3 being: 



_ 4K(m) ^^ 2 _cn(/3,m) 



N„ 



dn {(3, m) 



(22) 



Here, sn(a;, m), cn(a;, m), and dn(x, m) are the Jacobi 
elliptic functions of modulus m, while K(m) is the com- 
plete elliptic integral of the first kind |l6( . The two solu- 
tions have a common limit for m — > 1 giving the type III 
solution: 



sinh(/3) 



-^ t+5 \ (7V p ^oo), (23) 



cosh(/3(n + c)) 
with (3 being determined by: 

v = 2 cosh(/3) . 

A. Standard Poisson brackets, A = 



(24) 



In this case of the standard PB (A = 0), the equation of 
motion becomes: 

ilpn + bPn+l ~ 2V>™ + 1pn-l] 



v'(l+m n \ 2 )-vf'(\Tpn\ 2 )) 



Tpn = 0, (25) 



(l + /(l^l 2 )) 

or equivalently, 

(np n -2i> n + v'ip n )(l + f{\ip n \ 2 ) 

+ [V>„+1 + ^n-l] (1 + f(Hn\ 2 )) - Vf'(\i>n\ 2 )i>n = 0(26) 
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This equation has the stationary solutions of the form 
(Type I, II, and III) as given by Eqs. JEJ, and Q 
where the frequency lo is still given by Eq. (fTTJl but with 
v replaced by v'\ 



u> = 2 — v . 



(27) 



Also, the two equations needed to determine m and (3 are 
still the same. The condition on the amplitude is that 
A = \Jp + 1 and that the coefficients of the polynomial 
a.j in Eq. I|15|l are given by 



P 



n£>-fc) 



1 < 3 <P- 



3 (j + iy. (p+iy 

This specifies the equation completely. We find: 



1 + f(x) = 1 + 



P+i 



p+1 

For the first few values of p one finds 



3 , 27 1 



e x for p 



(28) 



(29) 



V = 


1 : 


OL\ 


p = 


2 : 


a i 


p = 


3 : 


ot\ 


p = 


4 : 


a i 



J) a 2 lpi a 3 2^6' 



a 2 



2', 



125 : 



tt4 



1 

3125 ' 



We note that p = gives the results of Ref. We 
also note that the p > cases are just a rescaling of the 
results of Ref. [l2( . Thus the higher order nonlinearities 
do not give anything essentially new in the A = case. 



B. Non-standard Poisson brackets, A = 1 

This case of the non-standard PB represents a general- 
ization of Eq. 0. From Eq. (0J and Eq. lfTS|) we get 

itp n + [lp n+1 + l/'n-l] [1 + f(Hn\ 2 )} 
+ (*/ - 2)[1 + f(\tPn\ 2 Mn ~ Vf'{\i>n\ 2 )i>n = . (30) 

In order to use the identities for the Jacobi elliptic func- 
tions we must choose v' = 2, and we get a restricted 
generalization of Eq. Jfjjl 



iijj n + [ip n+1 + l/'n-l] [1 + f(\^n\ 2 )} 

-Vf'{\^n\ 2 )^n =0. 



(31) 



This equation has the stationary solutions of the form 
(Type I, II, and III) as given by Eqs. JHJ, |2J, and O 
except here the equations determining m and (3, are 







2K{m) 



V — LO 



^ dn(/3, m) 
cn 2 (/3, m) 



(32) 



for the type I solutions, 



4K(m) cn(/3, m) 

f.J = 77 , V — LO = 2- 



N p ' " dn 2 (/3,m) 
for type II, and finally 

v — lo = 2 cosh(/3) 



(33) 



(34) 



for type III. The condition on the amplitude is now that 



A = ^Jp + iy tt^-j and that the coefficients of the poly- 
nomial ctj are given by 



.,= , 1 ,. n ^- fc V— i if 



J (7 + 1)! (p+1)^ 1 

i<i<P- 



(35) 



Note that the frequency lo now also appears as a param- 
eter. This completely specifies our equation. 
We find: 



,. /•„„ If" 1 "- n iHi ,y " ■'" 



f P+1 



p+1 



V — LO p 

In addition, for p — > oo 



w p+1/ ^ — a; x 
i 1 1 



v p+1 



! + /(*) 



LO 



and 



/'(aO-e-S-*. 
For the first few values of p we find: 



(36) 



(37) 



(38) 



p = 


1 : 


«i 


p = 


2 : 




p = 


3 : 


ai 









V — LJ ( I 1/ — LO \ 

2u \ 2v h 

V — LO /-j V — LO \ 

2v V 3f /' 

^ — LO (A V — LO \ 

2v \ 4f II 
(i/—io')3 /-i a; \ 

192i/3 V 4i/ /■ 



fl'2 
fl'2 



_ {v-lo)2 /-, 



3i/ I' 
j—wyz r-t _ v—u \ 

\2v2 \ X Au II 



Note that if lo = (i.e., = 1) then, as expected, 

ctj are the same in non-standard Poisson bracket case as 
in the standard Poisson bracket case. We would like to 
remind that if one also chooses v = 2 then one obtains 
the AL model and its higher order generalizations. 
It may be noted here that in order to obtain the above so- 
lutions we have made use of two identities for the Jacobi 
elliptic functions [l7|. The first identity is 

dn 2 (x)[dn(x + a) + dn(x — a)] — Adn(x) 
+.B[dn(a; + a) + dn(x - a)] , (39) 



where 



+ = 2ns(a)ds(a) , B = -cs 2 (a) 



(40) 



5 



Here ns(o) = 4^' ds ( a ) = SS' cs ( a ) = We 
have suppressed the modulus m in our notation here. 
On repeatedly multiplying both sides of the identity i|39|) 
by dn 2 (x) and simplifying yields the following general 
identity of arbitrary (odd) rank: 

dn 2n (x)[dn(a; + a) + dn{x - a)} = B n [dxi(x + a) 



(*)> 



(41) 



+dn(.T -a)]+Aj2 B^^dn 2 ^-^ 1 
3=1 

which has been used in deriving some of the above solu- 
tions. Here A, B are the same as given by Eq. (I40|l . 
The second identity we have used above is |17j 

mcn 2 (x)[cn(a; + a) + cn(a; — a)] = Acn(x) 
+B[cn(x + a) + cn(x - a)} , (42) 

where 

A = 2ns(a)cs(a) , B = -ds 2 (a) . (43) 

On repeatedly multiplying both sides of this identity 
by mcn 2 (x) and simplifying yields the following general 
identity of arbitrary (odd) rank: 

m™cn 2 "(a;)[cn(a; + a) + cn(a; - a)] = B n [cn(x + a) 

n 

+cn(x - a)} + AJ2 m n - j B : >- 1 cu 2{n -^ +1 (x) , (44) 
i=i 

which has been used in deriving some of the above solu- 
tions. Here A, B are the same as given by Eq. (|43[1 . 



C. The p — > oo limit and the linear limit 

Both A = and A = 1 cases have the same linear (small 
signal) limit: 

iip n -(2 + v- v')^ n + (^„+i + Vn-i) = 0. (45) 

For the A = case, the p — > oo limit gives the same 
equation f4*5|l . because in this limit l+/(a;) = f'{x) = e x . 
For the A = 1 case, however, 1 + f(x) and f'(x) differ 
as is clear from Eqs. (I37JI and (|38J) . Inserting these two 
equations into Eq. (|31|) yields 



#71 



lp n+ l + tp n -i 



+ ve^ { ^ n+ l + _l n - X ~ ih) ^ 0. (46) 

This equation is a little bit tricky since it contains terms 
of the order ^/p and v / pe p . We must balance the ^/pe p 
terms first giving: 

1pn+l + 1pn-l = {v - U)lp n , (47) 

and next balance the smaller terms giving 

iipn = coip n . (48) 

Note that in this case we end up with two coupled linear 
equations, one describing the spatial variation and one 
describing the temporal variation. 



IV. PROPERTIES OF THE NEW SOLUTIONS 
A. PN Barrier 

In a discrete lattice there is an energy cost associ- 
ated with moving localized modes (such as a soliton or 
a breather) by a half lattice constant. This is the cele- 
brated Peierls-Nabarro (PN) barrier [TsL IT^ | . As is well 
known, while for the AL (i.e. p = 0, A = 1, v = 2) case, 
this barrier is known to be zero, as shown by us in Ref. 
[12j, this barrier is nonzero in the case of the saturated 
DNLS model (i.e. p = 0, A = 0). It is then of significant 
interest to know whether this barrier exists in models 
with higher order nonlinearities. Since for p = 0, A = 0, 
we have already studied the various properties of the so- 
lutions in Ref. |l2j | and since higher order nonlinearities 
do not give anything essentially new, we will only study 
the properties of the A = 1 solutions. In particular, we 
will show that like the AL case, even for all the higher 
order models, the PN barrier is zero. 

In view of v 1 = 2, the Hamiltonian (|15|l (for A = 1) 
takes a simple form 



JY 



n+l 



C^+i + ^ln(l + /(|^n| 2 )) 



H 1 +H 2 , (49) 



while the power P is as given by Eq. HEJ). The first part 
of the Hamiltonian without the logarithmic term ( Hi ) is 
easily evaluated in the case of all the three solutions as 
given by Eqs. 1|19|) . (|21|l and (|23|l and it is easily shown 
that for all the three solutions the answer is independent 
of the constant c, i.e., the distance from a lattice point 
where the center of the elliptic soliton solution is located. 
For example, for solution of type I, we obtain 



H 



2A 2 N, 



1 cs 2 (f3,m 



^-[dn(/3,m)-cs(/3,m)Z(/3,m)], (50) 



while for the solution of type II 



2A 2 N 

H i = j--;— ^r-[mcn(/3, m)-ds(/3, m)Z{fi, m)\ , (51) 

cs z (p, m) 

and for solution of type III 

H[ u = -4 J 4 2 sinh(/3). (52) 

Here Z((3,m) is the Jacobi zeta function While de- 
riving these relations, use has been made of the identities 
in Ref. [13 

dn(x)dn(x+a) = dn(a)— cs(a)Z(a)+cs(a)[Z(x+a)— Z(x)] 

(53) 

mcn(x)cri(x + a) = ?7jcn(a) — ds(a)Z(a) 
+ds(a) [Z(x + a) - Z[x)} , (54) 
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We are unable to evaluate the expression for power 
P and (hence) the second term in the Hamiltonian H 2 
analytically for any nonzero p in the case of solutions of 
type I and II. However, this is easily accomplished in the 
case of the localized solutions (of type III). However, even 
without the explicit computation, it is easy to show that 
the PN barrier is zero in the case of all three solutions. 
Let us first explain the key idea. The power P is given 
by 



Unfortunately for the solutions of type I and II we are 
unable to perform the additions analytically and hence 
compute P and H 2 analytically. However, for the spa- 
tially localized solutions of type III this is easily accom- 
plished. We observe that for the localized solutions of 
type III, each term in the sum has the form 



ln[l + asech 2 /3(n + c) 



(61) 



JV 



p = J2 h ^ + € + ^i + 



o+2l 



^pYn 



(55) 



and as is well known from the AL case |2jj, this sum is 
c-independent and given by 



where ip n = </>„e~ l (" t+l5 ) and </>„ is easily obtained from 
Eqs. (HHJ), (211 an d ©■ The key point to note is that the 
expression under the logarithm can always be factorized 
as 



1 



1 



2ct\v 



(p + l)u 

which can be further factorized as 



{p+l)a p v 2p 
(56) 



1+ (p+l)iA ) [ 1 + a M 1 + a ^n]--^+apK], (57) 



where 01, a p are the roots of the above equation. 

For example 



"•• ^r~- 11 ■ 



(/ ' ' U "> M . (58) 



Hence the expression for power takes the simple form 

JV 



P = £ln 



n=l 
p N 



1 



(p+lK 



^^ln[l + a^ 2 J 
3=1 n=l 



(59) 



We now observe that in the celebrated AL case, the 
power P is given by 



JV 

E 

71=1 



ln[l 



(60) 



where <p n is either adn[/3(n + c — vt, to)] or dn replaced 
by cn[/?(n + c — vt, to)] or by sech[/3(n + c — vt)], and it is 
well known in that case that this sum is independent of c 
since P is a constant of motion and t and c always come 
together in this expression [2Cj . As a result, it immedi- 
ately follows that even for the higher order DNLS models 
the power P and hence H 2 must also be c-independent 
being a sum of p terms of the form same as that appears 
in the AL case. Thus we see that remarkably enough, 
even in higher order DNLS models the PN barrier is zero 
in the case of all three solutions. 



In [l + asech 2 f3(n + c)] = - [sum -1 y/K\ 2 . (62) 

P 



Thus, in principle we know P and H 2 for type III solu- 
tions. As an illustration, for the p — 1 case, it is easily 
shown that 



H 2 = vP = 2f3v 



2v 

+7 



sinh 



1v , v — uj, , „ 
1 sinh/3 



1 2 



(63) 



Generalizing, for arbitrary p, the power P and H 2 are 
given by 



H 2 = uP + 2(3v 
H — - > smh 



(p+i)- 



v — UJ 



-a,j sinh /3 



(.64) 



It is indeed remarkable that the PN barrier is not only 
zero for the AL model but for even higher order general- 
izations. It would be worthwhile to examine if the higher 
order models are also integrable, although perhaps the 
answer is likely not in the affirmative. 



B. Stability 

In order to study the linear stability of the exact solutions 
■0^ (j is I, II, or III) we introduce the following expansion 



1pn(t) =1p 3 n + Slp n {t)e 



(65) 



applied in a frame rotating with frequency u> of the so- 
lution. The stability analysis for the p = 0, A = case 
was carried out in Ref. [T^l and as seen above, higher or- 
der nonlinearities do not give any new solution. We will 
therefore only consider the stability of A = 1 solutions. 
Upon using this expansion into the equation of motion 
(for A = 1, v' = 2) 



itp n + [lp n+1 + 1p n -l][l + f(\^n\ 2 )] 



vf(\TPn\ 2 )*P n = 0. 



(66) 
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Next, retaining only terms linear in the perturbation, and 
taking into account the basic frequency u) of the unper- 
turbed solutions and the perturbations, we get 

i5i> n + (6ip n+ i + 8i/> n -i) [1 + /(IV^I 2 )] 

+ (W - vf{\^ n \ 2 ))8^ n + (r n ^n+l + Vn-O/'dVnl 2 ) 

-vf"{\^n\ 2 )\i> n \ 2 ) W„ + Sr n ) = 0. (67) 

Continuing by splitting the perturbation 5ip n into real 
parts Su n and imaginary parts 5v n {5ip n = 8u n + iSv n ) 
and introducing the two real vectors 



SU = {Su n } and SV = {Sv n } 

and the two real matrices A = {A nm } and B = 
by defining 

Anm = [6 n ,m+l + ^n,m-l][l + f(\^Pn\ 2 )} 

+ (u-vf{\^ n \ 2 )-2vf'{\i, n \ 2 )\^ 



., +V»-i]/'(|V'»| 2 ) 

l][l + /(l^l 2 )] 



(68) 



+ [W - v/'OVnl 2 )]^!, 



(69) 



where m ± 1 in the Kronecker <5 means: mil mod JV. 
Then Eq. I|67(l can be written compactly as 

-SV + A5U = 0, and 5U + BSV = 0, (70) 

where an overdot denotes time derivative. Combining 
these first order differential equations we get: 

SV + AB5V = 0, and 5U + BA5U = 0. (71) 

The two matrices A and B are symmetric and have real 
elements. However, since they do not commute AB and 
BA = (AB) T (T means transpose) are not symmetric. 
AB and BA have the same eigenvalues, but different 
eigenvectors. The eigenvectors for each of the two ma- 
trices need not be orthogonal. The eigenvalue spectrum 
{7} of the matrices AB and BA determines the stability 
of the exact solutions. If it contains negative eigenvalues, 
the solution is unstable. The eigenvalue spectrum always 
contains two eigenvalues which are zero. These eigenval- 
ues correspond to the translational invariance (c) and to 
the invariance of the solution tp n to a constant phase fac- 
tor e~ lS (i.e., translation in time), respectively. 

In Fig. 1 we show three examples of such stability eval- 
uation. Figure 1 shows the lowest (non-zero) eigenvalue 
from the spectra of type III solutions obtained for p = 
(solid line) , p — 1 (dashed line), and p = 2 (dashed- 
dotted line). For the integrable AL case (p = 0) we 
observe the expected result that all eigenvalues have a 
positive real value, indication that the exact solution is 
stable for all widths j3. In contrast, we see that for both 
p = 1 and p = 2 the solution becomes unstable for cer- 
tain values of (3. This instability occurs for relatively 
large values of (3, and thus when the solutions are very 
localized. 
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FIG. 1: Lowest non zero-eigenvalues for type III solutions 
[Eq. (23)] for the (solid line) integrable AL equation (p=0), 
and for the nonintegrable cases of p=l (dashed line), and p—2 
(dashed-dotted line). For p = 1 and p — 2, the eigenvalues 
are observed to become negative, indicating instability, for 
certain values of (3. 



DISCRETE BREATHERS AND MOVING 
SOLITONS 



In the following we focus on the AL limit of Eq. jST} 
and investigate the existence of localized stationary so- 
lutions, discrete breathers and moving excitations, by 
means of direct numerical integrations. In Fig. Ufa) we 
show the profiles of the localized states corresponding to 
the type III solution (soliton limit) for different values of 
p and for the same parameter f3 [notice that this parame- 
ter fixes the norm (power)]. We see that as p is increased 
the amplitude (and the norm) of the solution increases 
as a consequence of the higher order nonlincaritics. In 
panels (b), (c), and (d) of this figure we also show the 
time evolution of the stationary solutions for the cases 
p = 1,2,3, respectively (similar results are obtained for 
higher values of p). Notice that the states are quite sta- 
ble under time evolution, thus confirming the existence 
of stable localized solutions of the generalized Ablowitz- 
Ladik (GAL) equation with higher order nonlincarity. 

Next we concentrate on moving solutions and on dis- 
crete breathers. To this end, we recall that for the case 
p = (i.e., the usual AL equation) exact moving solutions 
of the traveling waves type are well known |14| . One can 
show that an extension of these p = moving solutions 
to the case p > is not possible if one assumes a traveling 
wave ansatz. The existence of a zero PN barrier, how- 
ever, strongly suggests the existence of moving solutions 
for all values of p. In the following we investigate this 
aspect by means of a direct numerical experiment. In 
this context, let us consider initial conditions which are 
a linear superposition of two exact stationary solutions 
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FIG. 2: Panel (a). Stationary solutions of the GAL equation 
(y — 2) for different values of p ranging from p — 1 (lower 
curve) top = 8 (upper curve). Other parameters are f3 = 0.25, 
m — 1. Panels (b), (c) and (d) show the time evolution of the 
solutions corresponding to the cases p — 1, 2, 3, respectively. 



of the form 

ip(n,t) = Ae~ iui 

dn[/3(n 



dn[/3(n - 
N 



N 



Xo)}- 



X )]e 



-iS 



(72) 



with A,lu given by our previous formulas. By properly 
choosing the initial distance 2Xq between the centers of 
the humps (in order to have a weak overlap) we can bring 
them in interaction and at the same time have a good 
initial condition which is very close to an exact solution. 
As for the usual NLS solitons, the interaction between 
the humps depends on their distance and on their phase 
difference S. In particular, we have that the two humps 
attract each other if they are in phase (6 — 0) and repel 
if they are out of phase (5 = tt). 

The existence of a zero PN barrier can then be checked 
by increasing the initial distance between two out-of- 
phase humps and observe if the humps are set in mo- 
tion by the mutual repulsion. Since by increasing the 
initial distance one considerably reduces the interaction 
between the humps (the interaction goes to zero expo- 
nentially with the distance) one has that for large sepa- 
rations motion can exist only if the PN barrier is zero. 
In Fig. |3]we show the results of such a numerical exper- 
iment by reporting the trajectories of the humps center 
(point of maximum amplitude) obtained from the GAL 
equation with p — 1 (AL with cubic-quintic nonlinear- 
ity) for different values of X$ and initial phase S = n. 
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FIG. 3: Time evolution of the position of the maxima of a 
two- hump solution of the GAL equation with p = 1, m = 1, 
originated from the initial condition in Eq. I72|l with the 
distance between the two humps increased in steps of 1 from 
Xo = 7 (lower slopes) to Xq = 10 (higher slopes). Other 
parameters are fixed as 8 — 7T, (3 = 1.25. 



We see that for small initial separations the two humps 
move in opposite directions with high velocity while as 
we increase the initial separation the velocity gets pro- 
gressively smaller. Our numerical investigations seem to 
indicate the absence of any critical threshold in the initial 
separation above which motion is stopped, which being in 
agreement with our analytical considerations about the 
absence of the PN barrier. 

This behavior is also seen from Fig. 0{a) in which the 
time evolution of two stationary solutions of the GAL for 
p = 1, initially displaced by a distance larger than their 
rest widths, is depicted. From this figure it is clear that 
there is practically no radiation generated during the mo- 
tion, thus making the hump dynamics very close to that 
of exact (traveling wave) solitons. In panel (b) of this 
figure we depict the time dependence of the amplitude 
during the hump motion of panel (a), from which we ob- 
serve that, except for the initial part where interaction 
dominates, a very regular pattern is generated. Notice 
that the amplitude oscillation lobes correspond to the 
vertical segments visible in the trajectories depicted in 
Fig. |21 the minima of the lobes corresponding to the 
times at which the maximum of the profile moves by one 
lattice site (i.e., to next vertical segment in the trajec- 
tory plot) . From this we infer that the reciprocal of the 
period of the oscillation in the hump amplitude is just 
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FIG. 4: Panel (a). Time evolution of two solitons of the GAL 
equation obtained for p = 1, m = 1 and initial condition 
taken as in Eq. 1721 with Xq = 4. Other parameters are 
fixed as S = it, f3 = 1.25. Panel (b). Time dependence of 
the humps amplitude (modulo square of the maximum of the 
profiles) depicted in the left panel. The minima correspond 
to the times at which the amplitude moves to the next lattice 
site. Parameters are the same as for panel (a). Panel (c). 
Time evolution of Eq. 1311 with p = 1, N = 100 and initial 
condition taken as in Eq. 1721 with Xo = 1, m = 1, /3 = 1.25 
and 5 = it. Panel (d). Same as in panel (c) but for p = 2 and 
X = 3. 



the hump velocity in lattice site units. It is remarkable 
that the absence of radiation in the hump dynamics is 
observed also for very long times and can survive multi- 
ple reflections. This is shown in panel (c) of Fig. 0]where 
the dynamics of two humps moving with a higher velocity 
(the velocity is increased by reducing the initial separa- 
tion) is depicted. Notice that the humps undergo several 
collisions with almost no radiation generated. This be- 
havior is a direct consequence of the zero PN barrier and 
is reminiscent of their soliton behavior. In panel (d) of 
Fig. 21 we show the same type of behavior for the case 
p = 2. This result indicates that moving localized states 
of the GAL equation may exist also for higher values of 
p (the stability of these solutions, however, may become 
critical as p is increased). The existence of a zero PN 
barrier makes it also possible to have discrete breathers 
in the GAL equation. In order to show this we repeat the 
same type of experiments considered above but with an 
attractive instead of a repulsive interaction, i.e. we take 
the initial humps to be in phase instead of out of phase. 

In Fig. [S] we depict the dynamics of two stationary 
solutions of the GAL equation for p = 1 and with a small 
value of (3 so that the initial profiles are very wide and 



FIG. 5: Left panel. Time evolution of Eq. (1311 with p = 1, 
N — 100 and initial condition taken as in Eq. I|72|l with 
Xq = 5, m = 1, P = 0.15 and 5 = 0. Right panel. Same as 
for the left panel but for p = 2 and Xo = 8. 



have a large overlap. We see that, due to the mutual 
attraction, the two localized solutions undergo breath- 
ing oscillations with apparently no radiation generated. 
Notice that the norm, which is fixed by the parameter 
(3, is below the instability threshold and the oscillations 
continue forever. This solution represents therefore an 
almost exact discrete breather of the GAL equation with 
higher order nonlinearity. In the right panel of Fig. [S| 
we show a discrete breather for the case p = 2, indicat- 
ing that these solutions may also exist with higher order 
nonlinearities. By increasing the norm or by increasing 
p, however, we find that an instability in the dynamics 
may develop at later times and the solution may become 
unstable. A preliminary investigation indicates the exis- 
tence of a critical threshold (which depends on p) below 
which stable stationary humps and discrete breathers are 
stable. This stability threshold may be reminiscent of 
the collapse threshold that exists in the continuous NLS 
with higher order nonlinearities, for norms (powers) ex- 
ceeding a critical value. For most physical applications, 
however, only the lowest higher order nonlinearities will 
be of interest (i.e., the cases p=l,2). In these cases sta- 
ble discrete soliton-like and discrete breather solutions 
are found to be stable for a wide range of parameters 
(which one should be able to check with linear stability 
analysis) . 

We remark that the absence of the PN barrier, the pres- 
ence of discrete breathers and the absence of emitted ra- 
diation during the hump dynamics could indicate a pos- 
sible complete integrability of the GAL equation. Al- 
though this cannot be concluded without further analy- 
sis, we remark that the vanishing of the PN barrier is a 
necessary but not a sufficient condition for integrability. 
In this context we remark that for the nonintegrable dis- 
crete </) 4 models |21| and for discrete sine- Gordon chains 
[22| it is possible to have a zero PN barrier and radia- 
tionless moving kinks for particular values of parameters. 
The GAL equation could possibly be another example in 
which this phenomenon may occur. 
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VI. CONCLUSIONS 

We have introduced a class of discrete nonlinear 
Schrodinger equations with arbitrarily high order non- 
linearities which include as particular cases the saturable 
discrete nonlinear Schrodinger equation and the AL equa- 
tion. We have obtained three different types of exact 
solutions (both spatially periodic in terms of Jacobi el- 
liptic functions and their limiting hyperbolic case) for 
these models as well as that of a higher order generaliza- 
tion of the saturable discrete nonlinear Schrodinger equa- 
tion and the Ablowitz-Ladik e quat ion. We then studied 
the Peierls-Nabarro barrier [l8Lll9| for these solutions in 
various models and found that it is zero indicating that 
the soliton-likc solutions move without experiencing any 
effect of the underlying discreteness, which is quite re- 
markable. We also studied the stability of these solutions 
under small perturbation and found that they are robust 
as well as stable. Finally, we investigated the collision 
of two hump solutions (both in-phase and out-of-phase 



cases) and found that they collide and move without any 
radiation. The out-of phase case indicates the formation 
of discrete breathers. These results are strongly sugges- 
tive of the integrability of the models introduced here, al- 
though we did not attempt to prove this rigorously. Our 
solutions and related properties are likely to be useful in 
many physical contexts including optical waveguides 0], 
Bose-Einstein condensates in optical lattices U , and 
nonlinear optics in the context of photonic crystals 7] . 
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